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ABSTRACT 

The Planck Early Release Compact Source Catalog (ERCSC) has offered the first op- 
portunity to accurately determine the luminosity function of dusty galaxies in the very 
local Universe (i.e. distances ;$ 100 Mpc), at several (sub-)millimetre wavelengths, us- 
ing blindly selected samples of low redshift sources, unaffected by cosmological evolu- 
tion. This project, however, requires careful consideration of a variety of issues includ- 
ing the choice of the appropriate flux density measurement, the separation of dusty 
galaxies from radio sources and from Galactic sources, the correction for the CO emis- 
sion, the effect of density inhomogeneities, and more. We present estimates of the local 
luminosity functions at 857 GHz (350 /im), 545 GHz (550 fim) and 353 GHz (850 /im) 
extending across the characteristic luminosity L*, and a preliminary estimate over a 
limited luminosity range at 217 GHz (1382 /im). At 850 /mi and for luminosities L £ 

our results agree with previous estimates, derived from the SCUBA Local Universe 
Galaxy Survey (SLUGS), but are higher than the latter at L ;$ L+. We also find good 

agreement with estimates at 350 and 500 /im based on preliminary Herschel survey 
data. 

Key words: galaxies: luminosity function - galaxies: photometry - galaxies: starburst 
- submillimetre: galaxies 



1 INTRODUCTION 
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Our knowledge of the (sub-)millimetre luminosity function 
of galaxies has substantially improved in the last few years. 
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Before the launch of the Herschel Space Observatory (Pil- 
bratt et al. 2010), no blind sub-mm surveys of sufficient area 
were available, and estimates of the local luminosity func- 
tions had to rely on follow-up of samples selected at other 
wavelengths. The most notable example is the local 850 /jm 
luminosity function derived from the SCUBA Local Universe 
Galaxy Survey (SLUGS; Dunne et al. 2000) that provided 
SCUBA photometry of 104 galaxies drawn from the IRAS 
Bright Galaxy Survey (Soifer et al. 1989). Better constraints, 
particularly on the faint end of the 850 /jm luminosity func- 
tion as well as estimates of the luminosity function at many 
wavelengths in the range 60 /mi-850 /im, were obtained by 
Serjeant & Harrison (2005) by modeling the spectral en- 
ergy distributions (SEDs) of all 15,411 IRAS PSCz galax- 
ies (Saunders et al. 2000). The SEDs were constrained by 
all available far-infrared and sub-mm colour-colour relations 
from the SLUGS and elsewhere. The Herschel surveys have 
allowed the first determinations of local luminosity functions 
at 250, 350 and 500 fim based on complete samples of sub- 
mm selected galaxies (Dye et al. 2010, Vaccari et al. 2010, 
Dunne et al. 2011). The results are in generally good agree- 
ment with the estimates by Serjeant & Harrison (2005). 

The Planck (Planck Collaboration I, 2011) surveys have 
allowed the construction of the first all-sky catalogs homoge- 
neously selected at several wavelengths from 350 /im to 1 cm, 
the Early Release Compact Source Catalog (ERCSC; Planck 
Collaboration VII, 2011). At sub-mm wavelengths the dom- 
inant population of extragalactic sources consists of nearby 
dusty galaxies, mostly at distances < 100 Mpc, all with 
spectroscopic redshift measurements. These are almost ideal 
characteristics for estimating local luminosity functions. In 
contrast, because of the much smaller covered areas, to get 
sufficient statistics with Herschel surveys it is necessary to 
select galaxies up to z = 0.1-0.2 where evolutionary effects 
may be non-negligible, requiring model-dependent correc- 
tions. In addition, for many relatively distant galaxies only 
photometric redshifts are available. 

On the other hand, very low redshift galaxy samples 
(z « 0.1), like those provided by Planck surveys, have 
their own drawbacks that must be dealt with carefully. First, 
proper motions can give large contributions to the measured 
redshifts, hence making them unreliable as distance indica- 
tors; redshift-independent distance indicators need then to 
be used as far as possible. Second, we live in the outskirts 
of the Virgo super-cluster implying that Planck samples 
contain strong density inhomogeneities, while the standard 
method to estimate the luminosity function [the 1/V ma x 
method (Schmidt 1968)] assumes a uniform distribution of 
galaxies. Although the effect of inhomogeneities is substan- 
tially mitigated by the (almost) full sky coverage, it must 
be taken into account. 

In this paper we exploit the Planck ERCSC to estimate 
the local luminosity functions of star-forming galaxies at 
217 GHz (1382 /mi), 353 GHz (850 /mi), 545 GHz (550 /im) 
and 857 GHz (350 /mi). Our analysis provides a useful 
complement to the recent paper by the Planck collabora- 
tion (Planck Collaboration VII, 2012), in which the sub- 
millimetre to millimetre number counts and spectral indices 
of dust-dominated galaxies, as well as of extragalactic radio 
sources, selected in the ERCSC are presented. 

The paper is organized as follows. In §[2]we describe the 
selection of the samples, their completeness and the correc- 



tion of the Planck fluxes for CO-line emission. In §[3]we deal 
with the methodology used to measure the luminosity func- 
tion and present the results, that, in §[3] are compared with 
earlier estimates. Our main conclusions are summarized in 
§0 Throughout this paper we adopt a flat cosmology with 
fio.ra = 0.27 and H Q = 70kms~ 1 Mpc" 1 . 

2 THE SAMPLES 
2.1 The ERCSC 

The Planck ERCSC lists all high- reliability sources, both 
Galactic and extragalactic, based on mapping the entire sky 
once and 60% of the sky a second time in 9 frequency bands 
centered at 30, 44, 70, 100, 143, 217, 353, 545 and 857 GHz. 
At the frequencies of interest here, the Full Width at Half 
Maximum (FWHM) of the Planck beam is 4.68 arcmin at 
217 GHz (1382 /im), 4.43 arcmin at 353 GHz (850 /im), 3.80 
arcmin at 545 GHz (550 /mi) and 3.67 arcmin at 857 GHz 
(350 /mi) (Planck Collaboration I, 2011). Below 353 GHz the 
detected extragalactic sources are mostly radio loud Active 
Galactic Nuclei ( AGNs) , while at frequencies ^ 353 GHz 
they are mostly dusty galaxies (Planck Collaboration XIII 
2011; Planck Collaboration VII 2012). Since we are inter- 
ested in the latter objects we focus our analysis on the three 
highest frequency Planck channels although we could also 
derive an estimate, over a limited luminosity range, of the 
luminosity function of dusty galaxies at 217 GHz. 

The ERCSC offers, for each source, four different flux 
density measurements (Planck Collaboration VII 2011). 
One is that determined by the source detection algorithm 
(FLUXDET); the others are based on aperture photome- 
try (FLUX), on fitting the source with the Planck point 
spread function at the location of the source (PSFFLUX), 
and on fitting the source with an elliptical Gaussian model 
(GAUFLUX). As stated in the Explanatory Supplement to 
the Planck ERCSC (Planck Collaboration 2011) the ERCSC 
flux density values need to be corrected by factors depending 
on the spectral shape of the sources. We have adopted the 
multiplicative correction factors appropriate for a spectral 
index a = 3 (F v oc u a ), i.e. 0.896, 0.887, 0.903 and 0.965 
at 217, 353, 545 and 857 GHz, respectively, as given in the 
Explanatory Supplement. 

It is important to remember that, in building the 
ERCSC, the emphasis was more on source reliability than 
on photometric accuracy or completeness (Planck Collabo- 
ration VII 2011), and the absolute calibration of flux-density 
values was required to be accurate to within about 30%, al- 
though the signal-to-noise of the sources are much higher. To 
investigate the best Planck flux density measurements and 
their photometric quality we have exploited Herschel flux 
density measurements from the Herschel Reference Survey 
(HRS; Boselli et al. 2010), as given in Ciesla et al. (2012; 
149 galaxies at 857 GHz and 86 at 545 GHz), from the Key 
Insights on Nearby Galaxies Survey, KINGFISH (Dale et 
al. 2012, 33 galaxies at 857 GHz and 33 at 545 GHz), from 
the Herschel Virgo Cluster Survey (HeViCS; Davies et al. 
2012, 39 galaxies at 857 GHz and 19 at 545 GHz), and from 
the Herschel Astrophysical Terahertz Large Area Survey (H- 
ATLAS; Eales et al. 2010) as given in Herranz et al. (2012, 
11 galaxies at 857 GHz and 6 at 545 GHz). We have applied 
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Figure 1. Comparison of the 4 Planck flux density estimations at 545 GHz (550 (im) with Herschel measurements at 500 fim, colour 
corrected to 550 /xm using the spectral index determined from the Planck measurements for each individual galaxy at 545 and 857 GHz. 
The HRS flux densities (black) are from Ciesla et al. (2012), Kingfish flux densities (red) are from Dale et al. (2012), the HeViCS flux 
densities are from Davies et al. (2012) and the H-ATLAS flux densities (blue) are from Herranz et al. (2012). The 4 panels refer to the 
4 estimates provided in the ERCSC: 'FLUX', 'PSFFLUX', 'GAUFLUX' and 'FLUXDET' (from left to right), see text. The black solid 
lines correspond to a Herschel/ Planck flux density ratio of 1. The dotted curves are the linear least squares fits [eq. {TJ] for the full 
sample (coefficients A and B in Table [TJ and the vertical dashed lines correspond to the adopted 80% completeness limit (see § 12.31 1. 
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Figure 2. Comparison of the 4 Planck flux density estimations at 857 GHz (350 /im) with Herschel measurements at the same frequency. 
The symbols and the lines have the same meaning as in Fig. [JJ 



to the Herschel flux density measurements colour corrections 
by factors of 1.004 and 1.02 at 350 and 500 /im, respectively, 
appropriate for extended sources with dust temperature of 
20 K and dust emissivity index /3 = 1.5 (Table 2 of Ciesla et 
al. 2012). 

The Planck 857 GHz and Herschel 350 /im bands have 
(by design) very similar central wavelengths and the colour 
corrected flux densities in these bands can be compared di- 
rectly. The Planck 545 GHz (=550 /im) and Herschel 500 /im 
bands, on the other hand, are slightly offset. In order to com- 
pare the measurements in these two nearby bands a small 
colour correction is needed. In order to do this we extrap- 
olate the Herschel flux densities to 550 /im using the spec- 
tral index for each source as determined by the measured 



Planck flux densities at 857 and 545 GHz. For sources miss- 
ing Planck measurements at one of these 2 frequencies we 
have adopted the average value found for the others, i.e. 
a = 2.7. The results of these flux comparisons are shown in 
Figs.[T]and[2] 

In doing this exercise we need to take into account that 
Planck measurements may suffer from blending of close-by 
sources, individually detected by Herschel. Visually inspect- 
ing each HRS galaxy associated to a Planck source, Ciesla 
et al. (2012) found potential contamination of 11 sources 
out of the 155 in common with Planck at 350 /im. Two 
Planck sources in the HeViCS survey, out of 60 observed, 
are found to be blends of galaxy pairs well resolved by Her- 
schel: NGC4298 + NGC4302 and NGC 4567/8 (Davies et 
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al. 2012). One additional blend of two comparably bright 
galaxies (NGC 3719 + NGC 3720) was found by Herranz et 
al. (2012), out of 12 detected galaxies in the H- ATLAS equa- 
torial fields. In these cases, the Planck source has been split 
into its components and we have adopted the Herschel flux 
densities of the individual galaxies at 350 and 500 (im (the 
latter rescaled to 550 (tin). At longer wavelengths we have as- 
signed to the two components flux density values whose sum 
equals the measured Planck flux density at that wavelength 
and whose ratio is equal to the one measured by Herschel at 
500 (tin. A visual inspection of optical images of the ERCSC 
sources without Herschel observations showed that only few 
per cent of them are associated to galaxy pairs. Therefore 
we conclude that source blending is not a major problem for 
estimates of the luminosity functions using Planck data. 
We have performed a linear fit to the data in Figs.[T]and 

m 



Fu 



= A + B X Fpi a 



(1) 



The derived coefficients are listed in columns 2 and 3 of 
Table[T] As the ERCSC is expected to suffer from the 
Eddington (1913) bias at low flux densities, we have re- 
peated the calculation of the linear fits restricting our- 
selves to the sub-samples of galaxies with fpianck > 
3 Jy at 857 GHz (350 (im) , and with i*pi anc k > 2 Jy at 
545 GHz (550 /im); the corresponding coefficients, A* and 
B» , are also given in Table [1] The rms fractional differ- 
ences between Planck and Herschel flux densities, rms = 

(l/v^)[Etl(^Planck, i -F H e fS chel, i ) 2 /^la M k,i] 1/2 I *™ the 

two bright sub-samples are given in the last column of the 
Table. 

From this analysis we draw the the following conclu- 
sions: 

• As expected, all the 4 Planck measurements are system- 
atically higher than the Herschel measurements for flux den- 
sities below ;S 1.5 Jy (the best fit values of the coefficient A 

are all negatives). This is likely due to the Eddington (1913) 
bias. 

• At 857 GHz, FLUX is the Planck flux density measure- 
ment most consistent with Herschel data, as it provides the 
lowest dispersion around the Herschel measurements. How- 
ever, compared to GAUFLUX, it seems to slightly underes- 
timate the flux density of the brightest fluxes, corresponding 
to resolved nearby galaxies. This is probably due to the fail- 
ure of the adopted point source model. 

• At 545 GHz the Planck flux density estimates most 
consistent with Herschel results are GAUFLUX and 
FLUXDET, the former showing a slightly lower dispersion 
of Planck to Herschel flux density ratios for the bright sub- 
sample. Both FLUX and PSFFLUX underpredict the flux 
densities of the brightest (i.e. resolved) galaxies, with de- 
rived values of the slope B significantly higher than 1. 

In the light of these considerations, we took FLUX as our 
reference Planck flux density measurement at 857 GHz and 
used eq. {1} with the values of A and B taken from Table[T] 
(first and second columns) to bring it in statistical agree- 
ment with Herschel measurements before deriving the rest- 
frame luminosities. At 545 GHz we adopted the GAUFLUX 
values corrected again according to eq. fTJ. At 353 GHz 
and 217 GHz there are no P/ancfc-independent flux density 
measurements for a statistically significant number of local 
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Table 1. Coefficients of the linear least square relation between 
Planck and Herschel flux densities at 857 and 545 GHz for the 
whole sample shown in Figs[T] and [2] (A, B) and for the sub- 
samples with flux density greater than 3 Jy at 857 GHz and 
greater than 2 Jy at 545 GHz (A, and £?*). Also shown are the rms 
fractional differences between Planck and Herschel flux densities 
for the bright sub-samples (rms*). 



galaxies in the Planck sample (in fact, only a few galaxies 
in our 353 GHz Planck sample have SCUBA imaging data 
from SLUGS). Since the luminosity function is derived for a 
sub-sample of relatively bright galaxies (see § 2.3), and given 
that GAUFLUX seems to perform quite well at both 857 and 
545 GHz for the sub-mm brightest galaxies, we have decided 
to use that flux density estimate at 353 and 217 GHz. 



2.2 Selection of star- forming galaxies 

The first step towards the selection of dusty galaxy sam- 
ples is the identification and removal of Galactic and AGN- 
dominated sources. The contamination by Galactic sources 
was minimized by masking regions heavily affected by Galac- 
tic emissions. Planck Collaboration VII (2012) used still un- 
published Planck maps to cre ate the mask. Here we rely in- 
stead on the 100 (im map of ISchlegel. Finkbeiner. fc Davisl 
(|l998h that we smoothed with a beam of FWHM= 180'. Pix- 
els in the upper intensity quartile (in the unsmoothed map) 
are discarded as being contaminated by Galactic emission. 
The masked region includes 29.1% of the sky, thus leaving 
an area of 29,250 deg 2 for the measurement of the luminos- 
ity function. The numbers of Planck detections within that 
area are listed in Tabled 

In order to help identifying Galactic sources, the 
ERCSC includes two flags for each source: 'EXTENDED' 
and 'CIRRUS'. The flag EXTENDED is set to 1 when the 
square root of the product of the measured major and mi- 
nor axes of the source is 1.5 times greater than the square 
root of the product of the major and minor axes of the 
estimated Planck point spread function at the location of 
the source. Otherwise the flag is set to 0. Given the ar- 
cmin size of the Planck beam in the frequency range 857 
to 217 GHz (see §H3J, detections with EXTENDED=1 are 
likely to be associa ted with structure in the Galactic in- 
terstellar medium (|Herranz et al.l 120121 ). The numbers of 
Planck detections with EXTENDED=0 within the selected 
area are reported in Table[2] together with that of sources 
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Figure 3. Histogram of the flag CIRRUS values for the local galaxies in our sample with EXTENDED=0 (red histogram) and with 
EXTENDED=1 (blue histogram). The sub-sample of galaxies with EXTENDED=1 that were flagged as "contaminated by cirrus" is 
shown by the yellow histogram. 



with EXTENDED =1 within the same area. The flag 'CIR- 
RUS' quantifies how crowded the region around each Planck 
detection is, by providing the number of sources within a 2° 
radius from the source, in raw 857 GHz catalogues. The num- 
ber is normalized to a peak value of unity. The normalization 
factor is obtained from the number density of sources in the 
Large Magellanic Cloud region where the highest concen- 
tration of 857 GHz sources is observed. High values of this 
flag, e.g. CIRRUS > 0.125 (Herranz et al. 2012), may indi- 
cate that the detection has a Galactic origin or, if not, that 
the flux measurement is severely contaminated by Galactic 
dust emission. There are exceptions. In fact, some nearby 
galaxies have either EXTENDED =1 or CIRRUS>0.125 or 
both. One example is M81, which has EXTENDED^ 1 (at 
217 GHz) and CIRRUS=0.434. Planck Collaboration XVI 
(2011) found that about 20% of the objects in their sample of 
robustly identified nearby galaxies are classified as extended 
in the ERCSC. The most effective way to separate local 
galaxies from Galactic sources would be to visually inspect 
optical to infrared imaging data for all the Planck detec- 
tions. This is the kind of approach we have decided to adopt 
here. However, as this method is time consuming, we have 
carried it out on all the detections with EXTENDED=0 and 
only on a sub-sample with EXTENDED=1 (as only a small 
fraction of local galaxies are expected to be flagged as ex- 
tended in the ERCSC). We have decided to not exploit the 
flag CIRRUS for reasons explained in § 12.2.21 where we also 
describe the definition and the analysis of the sub-sample 
with EXTENDED= 1 . We first focus on the catalogue with 
EXTENDED=0. 

2.2.1 Sample with EXTENDED=0 

In order to pick up radio loud AGNs we have exploited the 
strong difference of their spectral shape compared with that 
of dusty galaxies at mm and sub-mm wavelengths. In this 
wavelength range the spectral index (a, with F v oc v a ) 
of radio-loud AGNs is generally £ while that of dusty 

galaxies is ^ 2.5. Since the former are increasingly impor- 
tant with decreasing frequency, we have first focused on the 
217 GHz sample and cross-correlated it with the ERCSC at 
100 GHz in order to estimate the 100-to-217 GHz spectral in- 
dex. We found 309 (out of 473) sources with F 1Q o/F 2 n > 1 
and dropped them from the sample as being dominated by 
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1382 


850 


550 


350 


"tot 


964 


1762 


3781 


6004 


«EXT=0 


473 


584 


946 


1887 


"EXT=1 


491 
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2835 


4117 


n EXT=0, local 


42 


220 


598 


1463 


"EXT=1, local 


11(9) 


33(20) 


55(21) 


182(56) 



Table 2. Total number of ERCSC sources, ntot, outside the 
adopted Galactic mask. Also shown are: the numbers of sources 
with EXTENDED=0 (n EXT =o) and EXTENDED=1 (n EXT=1 ), 
the number of identified galaxies at z < 0.1 with EXTENDED=0 
("EXT=o,local) and with EXTENDED=1 (n EXT=1J ocal)- In 
parenthesis are the numbers of the galaxies with EXTENDED=1 
and z ^ 0.1 that were kept for estimating the luminosity function 
(see text for details). 



radio emission. Note that this approach is different from 
that adopted by Planck Collaboration VII (2012), where the 
857-to-545 GHz and 545-to-353 GHz flux density ratios were 
used to separate radio sources from galaxies whose spectral 
energy distribution is dominated by thermal dust emission. 

As a counter-check, we cross-correlated our 217 GHz 
sample with the CRATES catalogue (Healey et al. 2007) 
adopting a search radius of 3'. This radius is somewhat 
larger, for these Planck channels, than the 1/2 FWHM of 
the Planck beam found by Planck Collaboration XIV (2011) 
to be sufficient to locate any related source. CRATES pro- 
vides a nearly uniform extragalactic (|&| > 10°) coverage 
for flat-spectrum (a > —0.5) radio sources brighter than 
65 mjy at 4.8 GHz. As expected, since radio sources detected 
by Planck are generally flat-spectrum, most of the sources 
dropped as synchrotron dominated have a bright CRATES 
counterpart. The cross-correlation with CRATES however 
misses a bunch of very bright steep-spectrum sources (e.g. 
3C 207, 3C 216, 3C 286, 3C 380) and includes some sources 
whose (sub)-mm emission is clearly dominated by dust (e.g. 
M82, M83, M104). On the whole, the cross-correlation with 
CRATES confirmed our previous conclusions and did not 
reveal any misclassified source. 

The remaining 473 — 309 = 164 sources were inspected 
individually. We used the interactive software sky atlas Al- 
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adin (Bonnarel et al. 2000) to inspect the available optical 
to far-infrared imaging data (e.g. IRAS maps) and to query 
the NASA/IPAC Extragalactic Database (NEDfl around 
the position of each source. We removed all those objects 
that did not have a nearby (z < 0.1) galaxy within the 
adopted search radius but, instead, satisfied one (or more) 
of the following conditions: (i) association with an extended 
clumpy /filamentary structure in the IRAS maps (if avail- 
able); (ii) presence of a bright (i.e. -Fi.4ghz > 0.1 Jy) radio 
source; (iii) presence of a group/cluster of galaxies (exam- 
ples are Abell 2218 and the Bullet cluster); (iv) the Planck 
sources is un-detected in IRAS maps. The last condition is 
meant to eventually identify and remove false Planck detec- 
tions that may have been induced by background fluctua- 
tions, which are expected to be significant in low-resolution 
surveys (see e.g. Negrello et al. 2004, 2005). We identified 122 
objects obeying at least one of the conditions listed above 
(we will refer to them as "contaminants"). Removing those 
sources from the catalogue left us with a final, cleaned, sam- 
ple of 42 local galaxies at 217 GHz with EXTENDED=0, all 
being either Messier objects or galaxies listed in the New 
General Catalogue (NGC). We then moved to the 353 GHz 
sample and cross-correlated it with the 217 GHz catalogue, 
using a 3' search radius. We removed all the "contaminants" 
in common with the 217 GHz sample and kept all the sources 
that were identified as local galaxies at that frequency. We 
then inspected each of the remaining unclassified Planck de- 
tections individually, in Aladin, and removed from the sam- 
ple all those objects satisfying at least one of the conditions 
previously illustrated. We followed the same approach at 
545 and at 857 GHz after cross-correlating the catalogues 
with the lower frequency samples in order to exploit all the 
information already in hand. 

The whole cleaning process, performed on the catalogue 
of Planck detections with EXTENDED=0, produced a sam- 
ple of 220 local galaxies at 353 GHz, 598 at 545 GHz, and 
1463 at 857 GHz (see Table[2]), the majority of which are 
either Messier or NGC objects. The visual inspection of 
multi-wavelength ancillary data for each individual object 
minimized the risk of misidentifications and automatically 
provided, through the NED, the information on the redshift- 
independent measurements of distance, needed to derive the 
local luminosity function. 

2.2.2 Sample with EXTENDED=1 

Table[2] shows that, particularly at 857 GHz and 545 GHz, 
the Planck detections with EXTENDED=1 are a factor 
> 2 more abundant than those with EXTENDED=0. Vi- 
sually inspecting all those sources would be a very ineffi- 
cient method to identify local galaxies, as the majority of 
the Planck detections with EXTENDED=1 are expected to 
be Galactic. Therefore we have first reduced the sample by 
performing an automatic cross-matching with the NED and 
by retaining for inspection only the Planck detections that 
happen to have at least one source with a measured redshift 
2 < 0.1 within a distance of 3'. A galaxy like the Milky 
Way, with a linear size of ~30 kpc, would be assigned EX- 
TENDED=1 for distances < 20Mpc or, equivalently, for 
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Figure 4. Postage stamp image of the galaxy UGCA-021 [z = 
0.006622), at 0.468 fim (left-hand panel) and at 100 (right- 
hand panel). This galaxy, indicated by the red circle, is flagged 
as extended in the ERCSC at 545 GHz and has CIRRUS=0.094. 
However the IRAS map shows that it lies behind a prominent 
structure in the Galactic interstellar medium. 
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Figure 5. Postage stamp image of the galaxy M82 (z = 0.00068), 
at 100 fim. The galaxy, indicated by the red circle, has EX- 
TENDED^ at 545 GHz and CIRRUS=0.5. The dashed red cir- 
cle marks the region of 2° radius centered on the galaxy. The 
high value of the flag CIRRUS is likely due to the presence, in 
the vicinity, of a prominent structure in the Galactic interstellar 
medium, seen in the top-right corner of the image. 



redshifts z £ 0.005. Therefore our choice of an upper limit 
of 0.1 in redshift is very conservative. The numbers of ob- 
jects left over after the cross-matching with the NED are 
reported in Table[2] Upon visual inspection of each individ- 
ual source in Aladin, to identify and remove radio sources 
and Galactic objects, we flagged as "contaminated by cir- 
rus" those galaxies that happened to lie behind a prominent 
filamentary structure in the IRAS 100 fim map. These galax- 
ies are likely to have their Planck flux density measurements 
severely affected by Galactic emission and therefore we have 
decided to drop them from the final sample. The number of 
galaxies with EXTENDED=1 actually used for deriving the 
luminosity function is shown in parenthesis in Table[2] One 
might have used the flag CIRRUS to identify local galaxies 
falling behind prominent Galactic structures, by imposing 
an upper limit on the CIRRUS values, for example 0.125, as 
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809 


2.39 


0.55 


217 


1382 


30 


197 


2.27 


0.44 



Table 3. Flux density limits (-Fn m ) corresponding to 80% com- 
pleteness and numbers of local dusty galaxies brighter than these 
limits, used for estimating the luminosity functions. The quanti- 
ties log F and <T[ og p are the best-fit values of the parameters that 
define the completeness function [eq. I|13|l ]. where the flux density 
is in mjy. 

suggested by Herranz et al. (2012). However Fig.[3] demon- 
strates that this approach is prone to misidentifications. In 
fact, the figure shows the histogram of the flag CIRRUS val- 
ues for the sample of local galaxies with EXTENDED=0 (in 
red) and that with EXTENDED=1 (in blue), while the yel- 
low histogram refers to the sub-sample of galaxies with EX- 
TENDED^ that were flagged as "contaminated by cirrus". 
While the latter objects are indeed those with the highest 
CIRRUS values at 217 GHz and 353 GHz, the situation be- 
comes less well defined at shorter wavelengths where several 
of them turn out to have CIRRUS ;$ 0.1. One example is 
shown in Fig. [4] On the other hand there are some nearby 
galaxies that, upon inspection of the IRAS 100 /im maps, do 
not seem to be contaminated by Galactic cirrus but, never- 
theless, have CIRRUS>0.2. Examples are M64, M82, M90, 
NGC7392. As illustrated in Fig.[5]for M82, the high CIRRUS 
value in these sources is probably due to the presence, in the 
vicinity of the galaxy (i.e. within 2°) but not directly on top 
of it, of a prominent structure in the Galactic interstellar 
medium. 



2.3 Number counts and completeness 

The completeness of our samples of local dusty galaxies can 
be tested by means of the differential source counts. Apart 
from the effect of inhomogeneities in the spatial distribu- 
tion of galaxies, the counts must have an Euclidean slope 
(i.e. dN/dS oc S~ 2 ' 5 ) and the onset of incompleteness is 
indicated by a decline below the Euclidean power law. The 
counts were estimated using a bootstrap resampling method 
that allows us to account for the uncertainties in the flux 
density measurements. In practise we have generated 1000 
simulated catalogues by resampling, with repetitions, the in- 
put catalogue. In each simulation we have assigned to each 
source a flux density value randomly generated from a Gaus- 
sian probability distribution with a mean equal to the mea- 
sured flux density and dispersion a equal to the associated 
error. The distribution of source flux densities derived from 
the simulations were binned into a histogram and the mean 
value in each bin was taken as the measurement of the num- 
ber count in that bin, once divided by the bin size and by the 
survey area. The errors on the number counts were derived 
assuming a Poisson statistic, according to the prescriptions 
of Gehrels (1986). 

The results are shown in Fig.[6] in the form of Eu- 
clidean normalized differential number counts and in Fig[7] 



as integral number counts. As expected, the counts exhibit 
an Euclidean slope at the brightest flux densities, followed 
by a downturn (or a flattening in the cumulative counts), 
which is due to the onset of incompleteness in the sample. 
The counts are compared with those estimated for dust- 
dominated galaxies by Planck Collaboration VII (2012), who 
used (i) a Galactic mask derived from Planck data, (ii) an 
automatic procedure to identify dusty galaxies based on 
Planck colours, (iii) completeness corrections obtained from 
extensive simulations, (iv) FLUX as the reference flux den- 
sity measurement at all frequencies. The agreement is gen- 
erally good, particularly at 857 and 545 GHz. At 353 GHz 
and at 217 GHz the faintest counts produced by the Planck 
collaboration lie above the extrapolation of our measure- 
ments. This may be due to the appearance of a strongly 
evolving population with a "super-Euclidean" slope. How- 
ever the possibility of either an overestimate of the correc- 
tion for incompleteness or of a contamination of the sample 
e.g. by Galactic and/or radio sources can hardly be ruled 
out at this stage. 

Fitting the brightest part of the Euclidean normalized 
counts (i.e. above 5, 2, 0.8, and 0.5 Jy at 857, 545, 353 
and 217 GHz, respectively; these flux density limits are in- 
dicated by the vertical dotted lines in Figs. |S] and with a 
straight horizontal line we get 486±32, 103±8, 14±1 and 
1.7±0.4Jy 1 ' 5 sr- 1 at 857, 545, 353 and 217 GHz, respec- 
tively. These values are fully consistent with those estimated 
by Planck Collaboration VII (2012) for dusty galaxies: 
627±152, 125±15, and 15±6 at 857, 545 and 353 GHz, re- 
spectively (at 217 GHz the Planck collaboration provides the 
Euclidean plateau level for synchrotron dominated galaxies 
only). 

The completeness as a function of the flux density was 
estimated as the ratio of the observed counts to those ex- 
pected from the best-fit Euclidean counts. It is well approx- 
imated by an error function 

C(F) = erf [(log F - log F) 2 /< F ], (2) 

where log-F and ai og p are free parameters whose best fit 
values for each Planck channel are given in Table [3] For 
the estimation of the luminosity functions we adopted the 
flux density limits corresponding to an 80% completeness. 
These limits are listed in Tableland represent observed flux 
densities, i.e. no correction for either CO-line emission (see 
next sub-section) or Eddington-bias [i.e. eq. ([TJ] was applied 
when measuring the number counts. In fact those corrections 
are only introduced when calculating the rest-frame lumi- 
nosities. In the same table we also give the corresponding 
numbers of local dusty galaxies used to estimate the lumi- 
nosity functions. All of them have a spectroscopic redshift. 
The completeness curves [eq. (| 13p ] are also used to compute 
the weight factor for each source, Wi = l/C(Fi), Fi being 
the flux of the i-th galaxy, to be used in the estimate of the 
luminosity function [see eq. © and eq. (113p j. 

By comparing our number counts with those predicted 
by Serjeant & Harrison (dashed lines in Figs. [6] and [7|), we 
find that the latter lie significantly below the former, partic- 
ularly at the longest wavelengths, implying that the Serjeant 
& Harrison local luminosity functions were also underesti- 
mated. This is not surprising since the Serjeant & Harrison 
estimates were obtained extrapolating the IRAS 60 fim data 
to 850 /im using the sub-millimeter/far- infrared colour re- 
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Figure 6. Euclidean normalized number counts of local star-forming galaxies at 857 GHz (350 ^m), 545 GHz (550 (im), 353 GHz 
(850 /im) and 217 GHz (1382 (im) (red dots with error bars). Also shown, for comparison, are the counts estimated by Planck Collabo- 
ration VII (2012; blue squares, corrected for incompleteness), as well as the predictions by Serjeant & Harrison (2005; dashed lines). The 
thick grey line is the result of the fit to the data points with flux density above the value indicated by the vertical dotted line, assuming 
that the differential counts have an Euclidean slope. 
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Figure 7. Integral number counts of local star-forming galaxies at 857 GHz (350 /im), 545 GHz (550 /im), 353 GHz (850 (an) and 217 GHz 
(1382 fim) (red dots with error bars). The dashed line is the prediction by Serjeant & Harrison (2005). The thick grey line represents the 
integral number counts derived from the fit to the bright tail of the differential counts shown in Fig. [6] 
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Table 4. Euclidean normalized number counts of local star-forming galaxies at 857 GHz (350 /im), 545 GHz (550 /im), 353 GHz 
(850 /^m) and 217 GHz (1382 (im), for flux densities greater than the values indicated by the vertical dotted lines in Figs. |61 and 
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lation derived from the SLUGS l|Dunne et al J 1200(f ). The 
latter is known to be biased against galaxies with large cold 
dust components (see Vlahakis et al. 2005 and the discussion 
in |4|, which are instead readily detected by Planck (Planck 
Collaboration XVI 2011). 

2.4 Correction of Planck fluxes for CO line 
emission 

All the Planck High Frequency Instrument (HFI) bands, 
except that at 143 GHz, include emission from the CO 
molecule in low redshift galaxies. This emission is a con- 
taminant for our purposes as we are concerned with the 
dust continuum emission from star forming galaxies. Since 
most of the sources in our samples lack CO data and we are 
dealing with the properties of the population as a whole, 
we apply to the Planck fluxes average correction factors for 
contamination from the various CO lines before deriving the 
corresponding rest-frame luminosities. 

The correction factors are based on the correlation 
between the total far-infrared luminosity (Lfir) and the 
CO line emission, Leo, and specifically on the corre- 
lation presented by Genzel et al. (2010) who found a 
close to linear relation between Leo and LFiR(50-300/im) 
for normal star-forming galaxies with an average ratio 
Lfir/Lc O(1 _ 0) = 27Z/Q/(Kkms _1 pc 2 ). The units of 
L'eo(i-o)> (Kkms _1 pc 2 ) can be converted to solar lumi- 
nosities with the relation given by Solomon et al (1997), 
Lco/Lq = 3.2 x 10- 11 (// rcst /GHz) 3 (L' co /(Kkms- 1 pc 2 ), 
to find Lfir/Lcou— o) — 5.55 x 10 5 with both luminosities 
in solar units. Although there is evidence that the correla- 
tion is different for merger systems with high values of Lfir 
(see e.g. Daddi et al. 2010) we did not feel that there were 
enough such systems in our sample to justify refinement of 
the corrections. 

In order to correct the Planck fluxes contaminated by 
CO(J = 2-1) (217 GHz band), CO(J = 3- 2) (353 GHz 
band), and CO(J = 5 - 4) (545 GHz bandfl, we assumed 
fixed line strength ratios. The CO{J = 2 - 1)/(J = 1 - 0) 
ratio in antenna temperature units, r%i, was taken as 0.8, 
which is a value typical of nearby galaxies (Leroy et al. 2009). 
The CO(J = 3-2)/(J = 1-0) ratio, r 3 i, was taken as 0.66 
based on the results of Yao et al. (2003). A similar median 
value, r3i ;lno dian — 0.7, was found by Mao et al. (2012; their 
mean value is however somewhat higher, r3i jmcan ~ 0.81), 
while Papadopoulos et al. (2011) find 7-31. mcan ~ 0.49. Fi- 
nally, r 51 = CO{J = 5 - 4)/(J = 1 - 0) was taken as 0.3 
(Papadopoulos et al. 2011). 

The procedure adopted was as follows. The IRAS 60 
and 100 /im fluxes (available for all our sources) were used 
to calculate the FIR (50-300/xm) luminosity, Lfir, in the 
same way as done by Genzel et al. (2010). This gives a pre- 
diction for the CO (J=l-0) luminosity. This was then scaled 
by the line ratios specified above, converted to flux units 
[Fco{j,j-\)/F GO (i-o) = rji(^/,j-i/^i_o) 2 ], to give the lu- 
minosity in the various other CO transitions. The measured 
Planck flux density was then used with the predicted CO 

2 No correction was made for the CO(J = 7 — 6) and (J = 
8 — 7) lines that contaminate the 857 GHz band as their estimated 
contribution is negligible. 



flux to calculate the equivalent width of the CO line. The 
ratio of this width and the bandwidth of the Planck band 
(taken to be 30% of the centre frequency, as indicated in the 
ERCSC Explanatory Supplement) was used as a measure of 
the fractional contribution of the CO line to the Planck flux. 
The average correction factors for CO contamination were 
0.91, 0.97 and 0.99 at, 217, 357 and 545 GHz, respectively. 
These values are consistent with the level of contamination 
estimated by Planck Collaboration XVI (2011). 

2.5 Distances and rest-frame luminosities 

Although a spectroscopic redshift is available for all the 
galaxies in our samples, this does not provide, in general, 
an accurate estimate of the source distance, as in the local 
Universe peculiar motions may introduce significant extra 
red- or blue-shifts in galaxy spectra. Redshift-independent 
distance estimates should therefore be used as far as possi- 
ble. Such estimates are available in the NED for the majority 
of the galaxies in our samples: 65% of them at 857 GHz, 88% 
at 545 GHz, 95% at 353 GHz and 98% at 217 GHz. For the 
remaining galaxies, distances were computed from the red- 
shifts, previously corrected to the reference frame defined by 
the Cosmic Microwave Background radiation using the NED 
on-line calculatoi__. We assigned to the redshift- dependent 
distances an arbitrary error of 30%. The distribution of the 
source distances at each frequency, shown in the top panels 
of Fig. [__ illustrates that these are truly local samples: the 
distances for the flux-limited sample used to estimate the 
luminosity functions are all ^ 100 Mpc (i.e. z ^ 0.02). 
The rest-frame luminosity of the i— th galaxy at the fre- 
quency of observation u h B is given by 

47rd 2 

Li = 77"7 777 — ( 3 ) 

(1 + Zi)k(Zi) 

where L; is the source flux at v = v b s , di is the source 
distance, Zi its redshift and fc(z;) = L(j/ b s (l + Zi)) / L(f b s ) 
is the k— correction. The latter is estimated using the SED 
of an Sd galaxy taken from the SWIRE template librarjQ 
(Polletta et al. 2007). Figure _J] shows the distribution of flux 
densities as a function of source luminosity. This is meant 
to illustrate the luminosity range covered by the samples of 
sources above the adopted flux density limit and how the 
luminosity bins are populated. The distribution of luminosi- 
ties as a function of the source distance is shown in the 
lower panels of Fig.__| In practice, the luminosity function 
is estimated from galaxies within a distance of 100 Mpc, as 
indicated by the vertical dashed line in Fig. [8] 



3 LUMINOSITY FUNCTIONS 

We have determined the luminosity functions using two 
different methods: the 1/V ma x method and the parametric 
maximum likelihood method. Below we provide a brief de- 
scription of the two approaches. The luminosity function 
is denoted as 4>{L) and is taken to represent the comoving 

3 http:/ /ned. ipac.caltech.edu/help/velc_help. html 

4 www.iasf-milano.inaf.it/~polletta/tcmplatcs/swire_templatcs 
.html 
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Figure 8. Top panels: distributions of source distances at 857 GHz (350 /im), 545 GHz (550 /tin), 353 GHz (850 /tm) and 217 GHz 
(1382 /mi) for the total samples of local star-forming galaxies (red histograms), for the sub-sample without redshift-indcpendent dis- 
tance measurements (yellow histograms) and for the sub-sample with completeness above 80 % and luminosity higher than the adopted 
minimum value (black line). The redshift corresponding to a given distance as determined by the cosmic Hubble flow alone (i.e. the 
one derived by neglecting peculiar motion, that we refer to as effective redshift, z e fj) is shown in the upper x-axis. Bottom panels: 
monochromatic luminosity as a function of distance at 857, 545, 353 and 217 GHz. The colour code is the same as in the top panels. 
The dashed line shows the luminosity that corresponds to the adopted flux density limit as a function of distance, while the horizontal 
dotted line indicates the adopted minimum luminosity. Galaxies in the samples used to compute the luminosity functions have distances 
< 100 Mpc (a limit indicated by the vertical dashed line). 
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Figure 9. Measured flux density as a function of the estimated monochromatic luminosity at 857 GHz (350 /mi), 545 GHz (550 /mi), 
353 GHz (850 /mi) and 217 GHz (1382 /tm). The horizontal dashed line indicates the flux limit above which the sample is 80% complete 
(see q [23t . 



number density of objects per unit logarithmic interval in 
luminosity, i.e. 



4>{L) 



dN 



dV dlogi ' 



(4) 



3.1 The 1/Vmax method 

The standard way of deriving the luminosity function is via 
the 1/V m ax estimator (Schmidt 1968; Avni & Bahcall 1980) 



= HLj) = 



N 3 



Q Alogi ^ V m 



(5) 



where £1 is the solid angle of the survey, Wi are the weight 
factors taking into account the incompleteness of the sam- 
ples (see § 12. 3p and the sum is over all the Nj sources with 
luminosity in the range [logij — AlogL/2, logi, + AlogL/2]. 
The quantity V m ax,» represents the (comoving) volume, per 
unit solid angle, enclosed by the maximum (comoving) dis- 
tance, r maXj i, at which the i-th object is detectable, given 
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Table 5. Local luminosity function at 857 GHz (350 /im), 545 GHz (550 /im), 353 GHz (850 Aim) and 217 GHz (1382 ^m) measured via 
the l/Vmax method, in units of Mpc - 3 dex . Luminosities are in units of W/Hz. 
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Figure 10. Local luminosity functions at 857 GHz (350 /mi), 545 GHz (550 /im), 353 GHz (850 /im) and 217 GHz (1382 /im) estimated 
with the 1/Vmax method. Also shown, for comparison, are the LF measured by Dunne et al. (2000) and Vlahakis et al. (2005) at 353 GHz 
(850 /im; blue squares and green squares, respectively), and by Vaccari et al. (2010) at 350 /an and at 500 /an (blue triangles). The 
Herschel luminosities have been converted from 500 to 550 /an by assuming a spectral index a = 2.7, the mean value found for local 
galaxies in the Planck sample. The dashed curves at 857, 545 and 353 GHz represent the predictions by Serjeant & Harrison (2005). 
The dashed curve at 217 GHz is the extrapolation of the Serjeant & Harrison prediction from 353 GHz assuming again a spectral index 
a = 2.7. The vertical dotted lines correspond to the adopted minimum and maximum luminosities (see text). 
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Figure 11. V/Vmex test at 857 GHz (350 /mi), 545 GHz (550 /im), 353 GHz (850 /tm) and 217 GHz (1382 /mi) for the flux-limited samples 
used to measure the luminosity function. The red dots, with error bars, are the mean V/V m ax values for various luminosity bins while 
the orange error bar, here arbitrarily placed at a luminosity of 10 26 W/Hz, represents the ±ltr uncertainty on the V/ V ma x value for the 
whole sample. The vertical dotted lines have the same meaning as in Fie. 1101 
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the survey flux limit, Fu m , i.e., for the Euclidean geometry 
that applies here 



Vmax.i — g r max,i) 

where 



1 + Zi 



(1 + Zi)k(zi)Li 



47rF lin 



1/2 



(6) 



(7) 



As stated before, when no redshift-independent distance 
measurements are available we used the distance computed 
from the redshift of the source, i.e. d; = cZi/Ho, where Zi has 
been corrected to the reference frame defined by the Cos- 
mic Microwave Background radiation. However, as shown 
in Fig. [5] above the adopted completeness limit almost all 
the sources have an estimate of the distance that is inde- 
pendent of redshift. In order to account for uncertainties in 
both flux densities and distances we adopted a bootstrap 
resampling method to derive the luminosity function. We 
have generated 1000 simulated catalogues by resampling, 
with repetitions, the input catalogue. In each simulation, a 
value of the flux density and of the distance is randomly 
assigned to each galaxy by assuming a Gaussian distribu- 
tion with a mean equal to the measured values and a equal 
to the quoted errors. The rest-frame luminosities, the com- 
pleteness correction factors and the maximum volumes are 
re-estimated from the simulated fluxes and distances. The 
luminosity function is then derived using eq. ([5]) consider- 
ing only those objects with simulated flux density brighter 
than the adopted flux limit. Finally, the mean value of the 
simulated luminosity functions in each bin is taken as the 
estimate of the luminosity function in that bin and the rms 
of the distribution is taken as the uncertainty. The lumi- 
nosity function estimated with the 1/Knax method is shown 
in Fig. [TU] and tabulated in TableJS] In a flux limited survey 
the faintest luminosities are undersampled and, as a conse- 
quence, the luminosity function derived from the bootstrap 
method exhibits an artificial turn-off at faint luminosities. 
Therefore we have adopted, in our analysis, a minimum lu- 
minosity, I/min, corresponding to that of a source located 
4 Mpc away from us and whose flux density is equal to the 
adopted flux density limit. This choice avoids the appear- 
ance of the turnoff. For similar reasons we excluded from the 
estimated luminosity function simulated luminosities higher 
than the maximum observed luminosity in the flux limited 
sample. The adopted minimum and maximum luminosities 
are indicated by the vertical dotted lines in Figs. llOH12l 

The 1 / V m ax estimator assumes a uniform spatial distri- 
bution of sources. In general this is not the case, and the 
effect of large scale structures may not average away, par- 
ticularly when the sampled volume is relatively small. In 
order to check if the distribution of the sources within each 
luminosity bin is consistent with uniformity, a V/Vmax test 
was performed. In practice for each object within a given lu- 
minosity bin we computed the ratio between the comoving 
volume up to the source position, Vi, and the corresponding 
maximum accessible volume V ma x,i- For a uniform distribu- 
tion (VyVmax) = 0.5 ± 1/V12JV, where N is the number of 
sources in the bin. In order to account for the uncertainties 
in both flux density and distance measurements, the V/ Vmax 
values have been computed for each of the simulated samples 
used to estimate the luminosity function and then averaged 



for each luminosity bin. The results are shown in Fig. 1111 
In each panel, the orange error bar, arbitrarily placed at 
a luminosity of 10 26 W/Hz, represents the ±lcr uncertainty 
on the (V/Vmax) value for the whole sample. The latter is 
consistent with 0.5 for all the four samples, although some 
small deviations from 0.5 are observed in some luminosity 
bins at 857 GHz. We conclude that, within 2a, the distri- 
bution of galaxies in our samples is consistent with being 
statistically uniform. In fact, as Planck covers most of the 
sky, local inhomogeneities are significantly diluted. 

3.2 Parametric Maximum Likelihood method 

A way to overcome the issue of large scale structure is to use 
maximum likelihood techniques, which can be either para- 
metric or non-parametric (Sandage, Tammann & Yahil 1979; 
Efstathiou, Ellis & Peterson 1988; Saunders et al. 1990). The 
main assumption behind these methods is that the relative 
densities of galaxies of different luminosities are everywhere 
the same, and the normalization of the luminosity function 
is a local measure of density. In practice this means that the 
observed luminosity function can be expressed as the product 
of a universal luminosity function, independent of position, 
and a position-dependent term that accounts for local inho- 
mogeneities (Yahil, Sandage & Tammann 1980; Yahil et al. 
1991). 

In this work we only implemented the parametric 
method (Sandage, Tammann & Yahil 1979, STY method), 
which has the advantage, compared with the non-parametric 
approach, that no binning of the data is required. How- 
ever an analytic form of the luminosity function must be 
assumed. We adopt a double power-law functional form: 



0(L|iV',L„a,/3) =JV' 



+ 



(8) 



as it provides a good fit to the local luminosity fu nction 
of IRAS-selected galaxies l|Serieant fc Harrisonl 12005:) . The 
double power-law has 4 free parameters in total: the nor- 
malization factor, iV', the faint end slope, a, the bright end 
slope, /3, and the characteristic luminosity, L*, at which the 
transition from the two slopes occurs. For a given set of 
values of the free parameters, the probability that a source 
is observed with rest-frame luminosity Li, and with an ef- 
fective redshift z c ft t i (i.e. the redshift that the source would 
have at the distance di if peculiar motions were negligible) 
is thus given by 



Pi = i 4>(Li\N', L+, a, P) -^r (z e ff,i) ^ 
N dzdi I 



(9) 



where dV c is the comoving volume element and L m i nj i is the 
minimum luminosity that can be observed at the distance 
and the redshift of the i— th source, given the flux limit of 
the sample, i.e. 



-£rnin,' 



47rd 2 Fi im 
(1 + Zi)k(zi 



(10) 



N is the total number of objects with L L m in,i and F ^ 
Fh m given the model LF 

/Zmax nr POO 

dz dzln( z) J d\o^{L\N',U, a ,p), (11) 

where A n i n = max[logL m i n ,i, \ogLn m (z)], and Ln m (2) is the 
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Figure 12. Local luminosity functions at 857 GHz (350 /im), 545 GHz (550 /an) and 353 GHz (850 /im) derived using the maximum 
likelihood parametric method (orange curve) compared with that measured via the l/Vmax method (red dots with error bars). The grey 
shaded regions correspond to the 68 per cent confidence interval for the parametric estimate. The predictions by Serjeant & Harrison 
(2005) are also shown for comparison (dashed curve). The meaning of the vertical dotted lines is the same as in Fig. 1101 
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Table 6. Parameters of the maximum likelihood luminosity functions [eq. (O] at 857 GHz (350 (im) , 545 GHz (550 (im) , 353 GHz (850 /im) . 
At 217 GHz (1382 (im) the very limited statistic does not allow us to derive any meaningful constraint on the model parameters. The 
quoted errors correspond to the 68% confidence intervals. 



luminosity corresponding to the flux density limit at redshift 
z, while 2 m i n and 2 max are the minimum and the maximum 
redshifts observed in the flux-/luminosity- limited sample. 
The probability defined by eq. ® is independent of the nor- 
malization TV' as the latter cancels out in the ratio between 
the differential and cumulative luminosity functions. This 
is what makes the maximum likelihood approach indepen- 
dent of local inhomogeneities under the assumption of a 
universal luminosity function. The probability of observing 
TV objects with luminosities {L\,L2, ...,Ljv}, with redshifts 
{zi, 22, 2jv} and with physical distances {di, da, djv} is 



TV 



TV! 



II' 



(12) 



where we have assumed that the number of detected sources 
follows a Poisson distribution with expectation value equal 
to TV defined by eq. (|11[) . In order to account for the incom- 
pleteness in our samples we follow Patel et al. (2012), by 
replacing the total number of objects (TV) with ^ Wi and 
introducing a weighting term, Wi/(w) ({w) being the average 
weight) as an exponent of the individual source likelihoods, 
such that 



PocTV 



n 



(13) 



By maximizing the probability P (i.e. the likelihood of the 
sample) one obtains the set of values of TV', L*, a and j3 
that make the observed data set the "most probable" one. 
Following the same bootstrap resampling method used to 



measure the luminosity function via the 1/V m ax method, we 
have estimated the maximum likelihood values of the free 
parameters for each of the simulated samples. The derived 
luminosity functions and their uncertainties are shown in 
Fig. 1121 They are fully consistent with those derived from 
the l/Vmax method. We have checked that the consistency 
with the 1 /VJnax results is mantained even if we use a power- 
law/Gaussian analytic models for the LF, like the one sug- 
gested by Saunders et al. (1990), 



0(L|TV',L*,a,<7) = TV 



'(e) 



exp 



2a 2 



■ log 



(' + £) 



(14) 



The agreement between the two estimates further supports 
the conclusion that inhomogeneities associated with large 
scale structure in the local Universe are not an issue in 
our analysis. Besides, overdensities associated with galaxy 
clusters in the local Universe are less pronounced in far- 
infrared/sub- mm surveys than in optical/near-infrared sur- 
veys. In fact such regions, especially their dense core, are 
dominated by ellipticals which have relatively low dust emis- 
sion. The maximum likelihood values of the parameters and 
their 68% confidence intervals, derived from the distribution 
of values produced by the bootstrap, are given in Table[6]for 
the 3 higher frequencies (857, 545 and 353 GHz). As illus- 
trated by Fig. [12] the luminosity range covered by the 217 



5 In this model the LF behaves as a power law for L < < and 
as a Gaussian in log L for L > > L* 
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Figure 13. Distribution of source (physical) distances in the 
SLUGS sample (blue histogram) and in the Planck flux-limited 
sample used to measure the luminosity function (red histogram) . 
The value of the redshift corresponding to a given distance {ef- 
fective redshift, see text) is shown in the upper x-axis. 

GHz sample is too small to allow a reliable estimate of the 
slopes a and /3. 

As a sanity check, we use eq. (6) of Planck Collabora- 
tion VII (2012) to derive the normalization of the Euclidean 
plateau of the number counts from the estimated luminos- 
ity function and compare it with the values quoted in $ 12.31 
The results, based on the best-fit parameters listed in Ta- 
bled are 509, 109 and 14 Jy 1 ' 5 sr" 1 at 857, 545 and 353 GHz, 
respectively, in excellent agreement with those derived from 
the observed number counts. 



4 COMPARISON WITH EARLIER 
ESTIMATES 

4.1 Comparison with Dunne et al. (2000) and 
Serjeant & Harrison (2005) 

Dunne et al. (2000) have carried out a follow-up program at 
850 /im with SCUBA - the SCUBA Local Universe Galaxy 
Survey (SLUGS hereafter) — of a sample of galaxies selected 
from the IRAS Bright Galaxy Sample (BGS; Soifer et al. 
1989). The BGS is complete to a 60 ^m flux limit of 5.24 Jy 
at |fc| > 30° and 8 > -30°. The SLUGS sample includes all 
BGS galaxies with declination in the range —10° < 5 < 50° 
and with velocities above w m j n = 1900 km s -1 (correspond- 
ing to a minimum redshift z m i n = u m m/c = 0.0063). The 
minimum velocity is chosen to exclude those galaxies whose 
angular size exceeds the SCUBA field-of-view (~2 arcmin- 
utes). In Fig.[T2] we compare the distribution of source 
distances in the SLUGS sample with that of the Planck 
flux-/luminosity-limited sample from which the luminosity 
function has been measured. The peak observed around 
17 Mpc is partially due to the Virgo cluster. However we 
have checked that the estimated luminosity functions do not 



change significantly if the region associated with the Virgo 
Cluster is masked off. 

The overlap between the Planck and the SLUGS sam- 
ples is poor (only 5 galaxies are in common) mainly be- 
cause of the low- z cut adopted by Dunne et al. (2000). 
The luminosity function derived by the latter authors is 
shown by the blue squares in Fig. 1101 It agrees with the 
Planck measurements for L353GHZ ~ 10 23 W/Hz but lies sig- 
nificantly below them at fainter luminosities. A likely ex- 
planation of the difference is the bias against cold dusty 
galaxies implicit in the SLUGS sample. In fact, the selec- 
tion at 60 fim combined with the imposed minimum red- 
shift tends to favor galaxies with relatively bright infrared 
luminosities which usually have warmer SEDs (e.g. Smith et 
al. 2011). This bias has been demonstrated by Vlahakis et 
al. (2005) and, more recently, by Planck Collaboration XVI 
(2011) who have performed a detailed analysis of the SED 
properties of a sample of low-redshift galaxies extracted from 
the ERCSC. We confirm those results by comparing the sub- 
mm/far-infrared colours of the SLUGS sample with those of 
the Planck sources used to derive the luminosity function 
at 850 fim. The results are shown in Fig.[2] together with 
the track of a grey-body with dust emissivity index /3 = 1.3 
and temperatures ranging from 20 to 50 K (in steps of 5K). 
The SLUGS sample has higher dust temperatures than the 
850 /J,m Planck-selected sample, and the difference is more 
pronounced for galaxies with lower 850 /im luminosities (yel- 
low dots). 

Serjeant & Harrison (2005) have gone one step further 
and derived the local luminosity function of galaxies at many 
wavelengths from 70 fim to 850 /im by modeling the SEDs 
of all 15411 galaxies from the redshift survey of the IRAS 
Point Source Catalogue (PSGz; Saunders et al. 2000). The 
PSCz survey covers 84 per cent of the sky to a depth of 
0.6 Jy at 60 /im and has a median redshift of 8400 km s" 1 (i.e. 
2 = 0.028). Starting from the IRAS measurements at 60 and 
100 fim, Serjeant & Harrison have exploited the strong cor- 
relation observed in the SLUGS between the far-infrared lu- 
minosity and the sub-mm/far- infrared colour to predict the 
sub- mm fluxes of each of the PSCz galaxies. As stated by the 
authors, these predictions are sufficient to define the sub-mm 
fluxes to within a factor of 2. The Serjeant & Harrison re- 
sults at Planck wavelengths are shown by the dashed curves 
in Figs. [6] and [7] for the number counts and in Figs. [TOl and 
I12l for the luminosity functions. Not surprisingly, since again 
we are dealing with an IRAS-selected sample, the agreement 
is good at high luminosities but PiancA;-based estimates are 
higher below the characteristic luminosity L*. 



4.2 Comparison with Vlahakis et al. (2005) 

To check whether the IRAS selection misses populations of 
sub-millimeter emitting galaxies Vlahakis et al. (2005) ob- 
served with SCUBA a sample of 81 galaxies selected from 
the Center for Astrophysics (CfA) optical redshift survey 
(Huchra et al. 1983). They found that the ratios of the mass 
of cold dust to the mass of warm dust for these galaxies is 
much higher than for the SLUGS galaxies, thus demonstrat- 
ing that the SLUGS was missing a significant population of 
galaxies characterized by large proportions of cold dust. In 
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Figure 14. Sub-mm/far-infrared colour-colour plane for the 
353 GHz (850 /im) selected Planck sources with completeness 
above 80% and with luminosity < 10 22 8 W/Hz (red dots) or 
> 10 22 ' 8 W/Hz (yellow dots). The stars represent the mean val- 
ues. The SLUGS sample (Dunne et al. 2000; blue triangles) is 
shown for comparison. The solid line is the track for a grey-body 
with dust emissivity index ft = 1.3 and temperatures ranging 
from 20 up to 50 K in steps of 0.5 K. 
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Figure 15. Sub-mm/far-infrared colour-colour plane for the 
353 GHz (850 /im) selected Planck sources with completeness 
above 80% (yellow dots) and for the sample of optically selected 
galaxie of Vlahakis et al. (2005; blue triangles). The solid line is 
the the same as in Fig. 1141 



fact, the 850 fim luminosity function derived by Vlahakis et 
al. is higher than the one derived from the SLUGS. 

The green dots in the 850 /im panel of Fig. 1101 show the 
luminosity function of Vlahakis et al. obtained with the 
method of Serjeant & Harrison, but using sub-mm/far-IR 
colour relations re-calibrated on the SLUGS and the op- 



tically selected samples togethei[f|. The result is in better 
agreement with the ERCSC luminosity function, and, as ex- 
pected, the sub-mm/far-infrared colours of the Vlahakis et 
al. sample are consistent with those of the 850 /im Planck- 
selected sample, as shown in Fig. 1151 



4.3 Comparison with Vaccari et al. (2010) 

Vaccari et al. (2010) have estimated the local luminosity 
functions of galaxies selected at 250, 350 and 500 ^m in the 
Lockman Hole and XFLS fields (14.7 deg 2 ) observed dur- 
ing the Herschel Science Demonstration Phase (SDP) of 
HerMES (the Herschel Multi-tiered Extragalactic Survey; 
Oliver et al. 2012). Their samples cover the redshift range 
< z < 0.2 and comprise 210 sources at 350 fim and 45 
sources at 500 ^m; at both wavelengths they have spectro- 
scopic redshifts for ~ 85% of their sources. 

We have converted to 545 GHz (550 /im) their 500 /im 
luminosity function by assuming a — 2.7, as this is the mean 
value of the sub-mm spectral index found from the analysis 
of the SED of local galaxies in the Planck sample (Clemens 
et al., in prep.). Their results, shown by the blue triangles in 
Fig. 1101 generally agree, within the uncertainties, with ours, 
even though Vaccari et al. did not take into account any 
evolutionary correction, while significant evolution between 
< z < 0.1 and 0.1 < z < 0.2 was reported, at 250 /jm, by 
Dye et al. (2010) and Dunne et al. (2011). 

The estimates of the local (sub-)mm luminosity func- 
tion presented here are the first derived from a complete 
all-sky (sub-)mm selected sample of galaxies. Together with 
the estimates of the luminosity functions of higher redshifts 
galaxies produced by Herschel (Eales et al. 2010b; Gruppi- 
oni et al. 2010; Lapi et al. 2011), they will provide a critical 
constraint for any model of galaxy evolution. 



5 CONCLUSIONS 

The Planck ERCSC has provided the first samples of 
truly local (distance ^ 100 Mpc) galaxies blindly selected at 

(sub)-mm wavelengths, large enough to allow us to obtain 
accurate determinations of the local luminosity function at 
857, 545 and 353 GHz down to luminosities one order of 
magnitude fainter than previous estimates and well below 
the characteristic luminosity L+. We have also obtained the 
first estimate, albeit over a limited luminosity range, of the 
local luminosity function of dusty galaxies at 217 GHz. To 
get there several steps were necessary: 

• Identification of the ERCSC flux density measurement 
most appropriate for our galaxies and application of suitable 
corrections, 

• separation of dusty sources from radio loud AGNs, 

• removal of Galactic sources, 

• determination of the completeness limits at each fre- 
quency, 



6 The 850 (im luminosity function derived by Vlahakis et al. using 
the method of Serjeant & Harrison is in perfect agreement with 
that directly measured from the optically selected sample alone. 
We decided to show the former because it has a better statistic. 
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• correction of Planck fluxes for the contribution of CO 
lines. 

Although spectroscopic redshifts are available for all galaxies 
in our samples, they may not be good distance indicators be- 
cause of the potentially important contributions of peculiar 
motions. Therefore we have used, as far as possible, redshift- 
independent distance estimates that are available for most 
sources. The effect of local inhomogeneities in the galaxy 
distribution was investigated by means of the U/U max test 
in different luminosity bins. We also used, in addition to 
the classical 1/Vmn* estimator, a parametric maximum like- 
lihood estimator of the luminosity function, that is immune 
to the effect of inhomogeneities. 

The derived luminosity functions are in good agree- 
ment, at high luminosities, with those obtained using follow- 
up observations of 60 /jm selected IRAS galaxies (Serjeant 
& Harrison 2005; Dunne et al. 2000). Below the character- 
istic luminosity, L t , however, our estimates are significantly 
higher, due to the bias against cold, generally low-luminosity 
galaxies, inherent in the SLUGS sample, and close the the 
estimate by Vlahakis et al. (2005) based on SCUBA follow- 
up of an optically selected sample. Our estimates are also 
consistent with those by Vaccari et al. (2010) based on data 
from the HerMES survey at 350 and 500 /im. 

A detailed analysis of the astrophysical properties (such 
as star-formation rates, dust masses, dust temperatures, in- 
frared luminosity function etc.) of the galaxies used here 
to derive the local sub-mm luminosity function will be pre- 
sented in a forthcoming paper (Clemens et al. in prep.). 
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